##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)
library(maftools)
library(copynumber)

##########################################################################################

option_list <- list(
    make_option(c("--sample_file"), type = "character") ,
    make_option(c("--seg_file"), type = "character") ,
    make_option(c("--class_order_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combine/"
    sample_file <- paste(work_dir,"/config/tumor_normal.class.list",sep="")
    seg_file <- paste(work_dir,"/titan/Titan_all_seg.final.tsv",sep="")
    class_order_file <- paste(work_dir,"/config/Class_order.list",sep="")
    ######
    images_path <- paste(work_dir,"/images/cnv_burden",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_file <- opt$sample_file
seg_file <- opt$seg_file
class_order_file <- opt$class_order_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

info <- data.frame(fread(sample_file))
seg <- data.frame(fread(seg_file))
class_order <- data.frame(fread(class_order_file))

t_logR <- 0.2

###########################################################################################
col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")
col <- col[1:4]

col_im <- brewer.pal(9,"YlGnBu")[6:8]
names(col_im) <- c("IM(IGC)" , "IM(DGC)" , "IM(IGC_DGC)")

###########################################################################################

dat <- seg
dat$Sample <- gsub( "S_" , "S" , dat$Sample )
info$Tumor <- gsub( "S_" , "S" , info$Tumor )
info$Normal <- gsub( "S_" , "S" , info$Normal )

dat$seg.mean <- dat$Median_logR
dat$loc.start <- dat$Start
dat$loc.end <- dat$End
dat$Tumor <- sapply( strsplit( dat$Sample , "_" ) , "[" , 1)
dat$Normal <- sapply( strsplit( dat$Sample , "_" ) , "[" , 2)

dat <- merge( dat , info[,c("Tumor" , "Class" , "Type")] , by = "Tumor" )

###########################################################################################
### 计算CNV改变负荷

dat$TCN <- 2*2^dat$seg.mean
dat$Length <- as.numeric(abs(dat$loc.start - dat$loc.end))

loss_region <- dat[which(dat$seg.mean <= -t_logR),]
gain_region <- dat[which(dat$seg.mean >= t_logR),]

result_loss <- c()
result_gain <- c()
for(Tumor in unique(dat$Tumor)){
  print(Tumor)
  class_s <- unique(dat[which(dat$Tumor==Tumor),'Class'])
  normal <- unique(dat[which(dat$Tumor==Tumor),'Normal'])
  type <- unique(dat[which(dat$Tumor==Tumor),'Type'])
  
  loss_rate <- sum(loss_region[which(loss_region$Tumor==Tumor),'Length'])/sum(dat[which(dat$Tumor==Tumor),'Length'])
  result_loss <- rbind(result_loss,data.frame(Tumor=Tumor,Normal=normal,Rate=loss_rate,Class=class_s,Type= type))

  gain_rate <- sum(gain_region[which(gain_region$Tumor==Tumor),'Length'])/sum(dat[which(dat$Tumor==Tumor),'Length'])
  result_gain <- rbind(result_gain,data.frame(Tumor=Tumor,Normal=normal,Rate=gain_rate,Class=class_s,Type= type))
}

## 确定上下阈值
min_burden <- -0.0001
max_burden <- 0.62

result_gain$Class <- factor(result_gain$Class,levels= unique(class_order$Class), ordered=TRUE)
result_loss$Class <- factor(result_loss$Class,levels= unique(class_order$Class), ordered=TRUE)

###########################################################################################
## 同一个人IGC和DGC的样本合并
## 分情况画
## IM + IGC
## IM + DGC
## 一个人的多个样本取均值
## loss
dat <- result_loss
dat_loss<- c() 
for(Normal in unique(dat$Normal)){
  print(Normal)
  tmp <- dat[which(dat$Normal==Normal),]

  tmp <- tmp %>%
  group_by( Normal , Class , Type ) %>%
  summarize( CNV_Burden_Loss = median(Rate) )

  tmp <- na.omit(data.frame(tmp))

  dat_loss <- rbind(dat_loss,tmp)
}

## gain
dat <- result_gain
dat_gain <- c() 
for(Normal in unique(dat$Normal)){
  print(Normal)
  tmp <- dat[which(dat$Normal==Normal),]

  tmp <- tmp %>%
  group_by( Normal , Class , Type ) %>%
  summarize( CNV_Burden_Gain = median(Rate) )

  tmp <- na.omit(data.frame(tmp))
  dat_gain <- rbind(dat_gain,tmp)
}

## IGC的负荷是否显著大于DGC的
# dat_gainGC <- subset(dat_gain,Class=="IGC" | Class=="DGC" )
#p <- wilcox.test(CNV_Burden_Gain ~ Class, data = dat_gainGC , paired = F , alternative = "less" ,exact=FALSE)$p.value

my_comparisons_1 <- list( 
  c(1, 2) , c(1, 3) ,
  c(2, 3)
)
## 扩增负荷
dat_gain$Class <- factor(dat_gain$Class,levels=c("IGC" , "DGC" , "IM"),ordered=T)

images_name <- paste(images_path,"/CNV_Gain_burden_all.pdf",sep="")
plot <- ggplot(data=dat_gain,mapping = aes(x=Class,y=CNV_Burden_Gain))+
geom_boxplot(lwd=1.5,aes(color=Class)) +
geom_jitter(position=position_jitter(0.2),aes(color=Class)) +
scale_color_manual(values=col) +
xlab(NULL) +
ylab('CNV Gain Burden')+
theme_bw() +
stat_compare_means(comparisons = my_comparisons_1) +
ylim(min_burden,max_burden) +
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
      legend.position ='none', # 隐藏图例
      legend.title = element_blank() ,
      panel.grid.major=element_line(colour=NA),
      legend.text = element_text(size = 8,color="black",face='bold'),
      axis.text.x = element_text(size = 10,color="black",face='bold'),
      axis.text.y = element_text(size = 10,color="black",face='bold'),
      axis.title.x = element_text(size = 10,color="black",face='bold'),
      axis.title.y = element_text(size = 10,color="black",face='bold'),
      axis.line = element_line(size = 0.5)) 
ggsave(file=images_name,plot=plot,width=4,height=6)

## 丢失的负荷
dat_loss$Class <- factor(dat_loss$Class,levels=c("IGC" , "DGC" , "IM"),ordered=T)

images_name <- paste(images_path,"/CNV_Loss_burden_all.pdf",sep="")
plot <- ggplot(data=dat_loss,mapping = aes(x=Class,y=CNV_Burden_Loss))+
geom_boxplot(lwd=1.5,aes(color=Class)) +
geom_jitter(position=position_jitter(0.2),aes(color=Class)) +
scale_color_manual(values=col) +
xlab(NULL) +
ylab('CNV Loss Burden')+
theme_bw() +
stat_compare_means(comparisons = my_comparisons_1) +
ylim(min_burden,max_burden) +
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
      legend.position ='none', # 隐藏图例
      legend.title = element_blank() ,
      panel.grid.major=element_line(colour=NA),
      legend.text = element_text(size = 8,color="black",face='bold'),
      axis.text.x = element_text(size = 10,color="black",face='bold'),
      axis.text.y = element_text(size = 10,color="black",face='bold'),
      axis.title.x = element_text(size = 10,color="black",face='bold'),
      axis.title.y = element_text(size = 10,color="black",face='bold'),
      axis.line = element_line(size = 0.5)) 
ggsave(file=images_name,plot=plot,width=4,height=6)

###########################################################################################
## 输出CNV的负荷的具体值
d1 <- dat_gain %>%
group_by(Class) %>%
summarize(Class=unique(Class), Sample_Num = length(CNV_Burden_Gain),
  Mean_Burden_Gain=round(mean(CNV_Burden_Gain),3),Median_Burden_Gain=round(median(CNV_Burden_Gain),3),
  Sd_Burden_Gain=round(sd(CNV_Burden_Gain),3))

d2 <- dat_loss %>%
group_by(Class) %>%
summarize(Class=unique(Class),Sample_Num = length(CNV_Burden_Loss),
  Mean_Burden_Loss=round(mean(CNV_Burden_Loss),3),Median_Burden_Loss=round(median(CNV_Burden_Loss),3),
  Sd_Burden_Loss=round(sd(CNV_Burden_Loss),3))

d3 <- merge(d1 , d2)
d3$Class <- factor(d3$Class,levels=class_order$Class)
d3 <- d3[order(d3$Class),]

images_name <- paste(images_path,"/CNV_Burden_ALL.tsv",sep="")
write.table(d3,images_name,sep="\t",quote=F,row.names=F)

###########################################################################################
## 不同亚型的肠化

my_comparisons_1 <- list( 
  c(1, 2) , c(1, 3) ,
  c(2, 3)
)

dat_plot_tmp <- subset( dat_gain , Class == "IM" )
dat_plot_tmp$Class <- ifelse( dat_plot_tmp$Type == "IM + IGC" , "IM(IGC)" , dat_plot_tmp$Class )
dat_plot_tmp$Class <- ifelse( dat_plot_tmp$Type == "IM + DGC" , "IM(DGC)" , dat_plot_tmp$Class )
dat_plot_tmp$Class <- ifelse( dat_plot_tmp$Type == "IM + IGC + DGC" , "IM(IGC_DGC)" , dat_plot_tmp$Class )
dat_plot_tmp$Class <- factor( dat_plot_tmp$Class , levels = c("IM(IGC)" , "IM(DGC)" , "IM(IGC_DGC)") , order = T )
images_name <- paste(images_path,"/CNV_Gain_burden_all.IM.pdf",sep="")
plot <- ggplot(data=dat_plot_tmp,mapping = aes(x=Class,y=CNV_Burden_Gain))+
geom_boxplot(lwd=1.5,aes(color=Class)) +
geom_jitter(position=position_jitter(0.2),aes(color=Class)) +
scale_color_manual(values=col_im) +
xlab(NULL) +
ylab('CNV Gain Burden')+
ylim( 0,0.03 ) +
theme_bw() +
stat_compare_means(comparisons = my_comparisons_1 , label.y = c(0.005,0.01,0.015)) +
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
      legend.position ='none', # 隐藏图例
      legend.title = element_blank() ,
      panel.grid.major=element_line(colour=NA),
      legend.text = element_text(size = 8,color="black",face='bold'),
      axis.text.x = element_text(size = 10,color="black",face='bold'),
      axis.text.y = element_text(size = 10,color="black",face='bold'),
      axis.title.x = element_text(size = 10,color="black",face='bold'),
      axis.title.y = element_text(size = 10,color="black",face='bold'),
      axis.line = element_line(size = 0.5)) 
ggsave(file=images_name,plot=plot,width=4,height=6)



dat_plot_tmp <- subset( dat_loss , Class == "IM" )
dat_plot_tmp$Class <- ifelse( dat_plot_tmp$Type == "IM + IGC" , "IM(IGC)" , dat_plot_tmp$Class )
dat_plot_tmp$Class <- ifelse( dat_plot_tmp$Type == "IM + DGC" , "IM(DGC)" , dat_plot_tmp$Class )
dat_plot_tmp$Class <- ifelse( dat_plot_tmp$Type == "IM + IGC + DGC" , "IM(IGC_DGC)" , dat_plot_tmp$Class )
dat_plot_tmp$Class <- factor( dat_plot_tmp$Class , levels = c("IM(IGC)" , "IM(DGC)" , "IM(IGC_DGC)") , order = T )
images_name <- paste(images_path,"/CNV_Loss_burden_all.IM.pdf",sep="")

plot <- ggplot(data=dat_plot_tmp,mapping = aes(x=Class,y=CNV_Burden_Loss))+
geom_boxplot(alpha =1 , size = 0.9 , width = 0.6 , outlier.shape = NA) +
geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
scale_color_manual(values=col_im) +
xlab(NULL) +
ylab('CNV Loss Burden')+
ylim( 0 , 0.01 ) +
theme_bw() +
stat_compare_means(comparisons = my_comparisons_1 , label.y = c(0.003,0.005,0.006)) +
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
      legend.position ='none', # 隐藏图例
      legend.title = element_blank() ,
      panel.grid.major=element_line(colour=NA),
      legend.text = element_text(size = 8,color="black",face='bold'),
      axis.text.x = element_text(size = 10,color="black",face='bold'),
      axis.text.y = element_text(size = 10,color="black",face='bold'),
      axis.title.x = element_text(size = 10,color="black",face='bold'),
      axis.title.y = element_text(size = 10,color="black",face='bold'),
      axis.line = element_line(size = 0.5)) 
ggsave(file=images_name,plot=plot,width=4,height=6)


###########################################################################################
## GC作为整体
my_comparisons_1 <- list( c(1, 2) )
## 扩增负荷
dat_gain$Class <- ifelse( dat_gain$Class != "IM" , "GC" , "IM" )
dat_gain$Class <- factor(dat_gain$Class,levels=c("IM" , "GC"),ordered=T)

images_name <- paste(images_path,"/CNV_Gain_burden_all.GC.pdf",sep="")
plot <- ggplot(data=dat_gain,mapping = aes(x=Class,y=CNV_Burden_Gain))+
geom_boxplot(lwd=1.5,aes(color=Class)) +
geom_jitter(position=position_jitter(0.2),aes(color=Class)) +
scale_color_manual(values=col) +
xlab(NULL) +
ylab('CNV Gain Burden')+
theme_bw() +
stat_compare_means(comparisons = my_comparisons_1) +
ylim(min_burden,max_burden) +
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
      legend.position ='none', # 隐藏图例
      legend.title = element_blank() ,
      panel.grid.major=element_line(colour=NA),
      legend.text = element_text(size = 8,color="black",face='bold'),
      axis.text.x = element_text(size = 10,color="black",face='bold'),
      axis.text.y = element_text(size = 10,color="black",face='bold'),
      axis.title.x = element_text(size = 10,color="black",face='bold'),
      axis.title.y = element_text(size = 10,color="black",face='bold'),
      axis.line = element_line(size = 0.5)) 
ggsave(file=images_name,plot=plot,width=4,height=6)

## 丢失的负荷
dat_loss$Class <- ifelse( dat_loss$Class != "IM" , "GC" , "IM" )
dat_loss$Class <- factor(dat_loss$Class,levels=c("IM" , "GC"),ordered=T)

images_name <- paste(images_path,"/CNV_Loss_burden_all.GC.pdf",sep="")
plot <- ggplot(data=dat_loss,mapping = aes(x=Class,y=CNV_Burden_Loss))+
geom_boxplot(lwd=1.5,aes(color=Class)) +
geom_jitter(position=position_jitter(0.2),aes(color=Class)) +
scale_color_manual(values=col) +
xlab(NULL) +
ylab('CNV Loss Burden')+
theme_bw() +
stat_compare_means(comparisons = my_comparisons_1) +
ylim(min_burden,max_burden) +
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
      legend.position ='none', # 隐藏图例
      legend.title = element_blank() ,
      panel.grid.major=element_line(colour=NA),
      legend.text = element_text(size = 8,color="black",face='bold'),
      axis.text.x = element_text(size = 10,color="black",face='bold'),
      axis.text.y = element_text(size = 10,color="black",face='bold'),
      axis.title.x = element_text(size = 10,color="black",face='bold'),
      axis.title.y = element_text(size = 10,color="black",face='bold'),
      axis.line = element_line(size = 0.5)) 
ggsave(file=images_name,plot=plot,width=4,height=6)
